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We study the critical behavior of the pair annihilation model (PAM) with diffusion in one, 
two and three dimensions, using the pair approximation (PA) and Monte Carlo simulation. 
Of principal interest is the dependence of the critical creation rate, Ac, on the diffusion 
probability in particular, whether survival is possible at arbitrarily small creation rates, 
for sufficiently rapid diffusion. Whilst the PA predicts that in any spatial dimension d > 1, 
Ac — at some diffusion probability D* < 1, Katori and Konno [Physica A 186, 578 (1992)] 
showed rigorously that for < 2, one has Ac > for any D < 1. Our simulation results 
are consistent with this theorem. In two dimensions, the extinction region becomes narrow 
as D approaches unity, following Ac oc exp[— const./(l — D)'*'], with 7 = 1.41(2). In three 
dimensions we find D* = 0.333(3). The simulation results confirm that the PAM belongs to 
the directed percolation universality class. 
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I. INTRODUCTION 



In recent decades, general theories of phase transitions and critical phenomena have 
been developed, unifying our understanding of equilibrium phase transitions in liquid- vapor, 
magnetic, liquid crystals and other systems. By contrast, the study of nonequilibrium 
critical phenomena is still in development. Since the transition rates in such systems do 
not satisfy detailed balance, the steady-state probability distribution in these systems is 
not known a priori, and the analysis must be based upon the dynamics. Starting from the 
basic contact process [l| , many particle systems have been studied in efforts to characterize 
scaling properties at nonequilibrium phase transitions {2 4|. These models, which involve 
creation and annihilation of particles on a lattice, typically exhibit a phase transition to 
an absorbing state (one allowing no escape), and so violate the detailed balance principle. 
An issue that has attracted some interest is the combined effect of multiparticle rules and 
diffusion (hopping), which tends to spread particles uniformly over the system. 

In this work we revisit the pair annihilation model (PAM) [s, 6|. In this model particles 
diffuse on a lattice at rate D, nearest-neighbor pairs of particles are annihilated at rate 
(1— D)/(l + A), and particles attempt to create new particles at rate {1—D)X/{1+X). Double 
occupancy of sites is forbidden. The model exhibits active and absorbing phases, separated 
by a continuous phase transition at Ac(-D). Using cluster approximations and Monte Carlo 
simulation, we determine the phase boundary in one, two, and three dimensions. 

The pair approximation predicts that for a diffusion rate greater than a certain value, 
D*, the critical parameter Ac = 0. (That is, for D > D*, an arbitrarily small creation rate 
is sufficient to maintain a nonzero particle density.) This prediction is known to be wrong 
in dimensions d < 2: Katori and Konno [3] proved that Ac > for any diffusion probability 
D < 1, in one and two dimensions. Their theorem is based on a relation between the PAM 
a.d*eb..a.chi.ga.„im.ati„g..a.doo,wa.kofB.a,„so.a„dGraya. Existence of a Z,' < 1 
is not ruled out in > 3 dimensions. The difference is connected with the nonrecurrence of 
random walks in d > 3 [7|. How Ac tends to zero as — ?■ 1 is, however, unknown. Moreover 
the question of whether, in three or more dimensions, D* is in fact less that unity, has not, to 
our knowledge, been studied. The principal motivation for the present work is to determine 
Ac(-D) via numerical simulation. We also verify that the model belongs to the directed 
percolation (DP) universality class, as expected on the basis of symmetry considerations 



[9|, llO| . Our simulation results, while consistent with the Katori-Konno theorem, show that 
in the two-dimensional case, Ac becomes extremely small as D approaches unity, possibly 
suggesting the (incorrect) impression that the critical value is actually zero at some finite 
diffusion rate. 

The remainder of this paper is organized as follows. In the following section (II) we 
define the model and discuss its limiting behaviors in the \-D plane. Then in Sec. Ill we 
present, for completeness, the one- and two-site cluster approximations. Simulation results 
are discussed in Sec. IV, followed by a brief discussion in Sec. V. 



II. THE MODEL 



The PAM is defined on a lattice, in which sites can be either occupied by a particle or 
vacant j^,!^,^]; we denote these states by cr^: = 1 (site x occupied) and cr^ = (site x vacant). 
There are three kinds of transition: nearest-neighbor (NN) hopping ("diffusion"), creation, 
and pairwise annihilation, with associated rates (1 — D)\/{1 + A), and {1 — D)/{1 + A), 
respectively. (Since the rates are parameterized so as to sum to unity, we are free to refer to 
D as the diffusion probability.) At each step of the evolution, the next attempted transition 
is taken as diffusion, creation, or annihilation, with probabilities D, (1 — D)X/{1 + A), and 
(1 — -D)/(l + A), respectively. 

In a hopping transition, a site x is chosen at random, along with a nearest-neighbor site 
y of X. Then if ax 7^ o-y, the states are exchanged. In a creation event, a site x is chosen. If 
ax = 1, a. nearest-neighbor y is chosen, and if ay = this variable is set to one. Finally, in 
an annihilation event, a pair of nearest-neighbor sites x and y are chosen at random, and if 
o'x = o'y = 1, both variables are set to zero. Each transition corresponds to a time increment 
At = 1/Nsite, where NsUe is the number of lattice sites. 

To improve efficiency, in simulations the site x is chosen from a list of occupied sites; 
then the time increment is At = 1/Np, with Np the number of particles in the system, 
immediately prior to the transition. In this implementation, the rate of annihilation of a 
given NN particle pair is 

~ At 1 + A iVp2^ ~dl + X 
where the factor 2/Np arises because either particle in the pair can be selected from the list 



of Np particles. 

The particle-free configuration is absorbing. By analogy with the contact process 
we expect that in the infinite-size limit the system undergoes a phase transition between an 
active state (with nonzero stationary particle density) and an absorbing one, as one crosses 
the critical line Ac(-D) in the X-D plane. As creation depends upon a single particle, the 
order parameter is the stationary density of particles, p. 

When a new particle is created, it always forms a pair with the "parent" particle, making 
these two particles susceptible to annihilation. In the active stationary state, increasing D 
at fixed A tends to reduce the fraction of nearest-neighbor particle pairs toward its random 
mixing value, p^. Thus we should expect Ac to be a decreasing (or at least, nonincreasing) 
function of D. In the simplest mean-field theory analysis, the annihilation rate is propor- 
tional to p^, so that for small p, one has p oc Xp — const, x p^, which admits a stationary 
solution p oc A for any nonzero creation rate. In the limit D — t- 1 we expect mean-field 
theory to hold, so that Ac — )■ in this limit. This raises the question of whether Ac vanishes 
at some diffusion probability D* strictly less than unity. While the two-site approximation 
predicts D* < 1 in any dimension, the results of Katori and Konno I?] imply that D* = 1 
in dimensions d <2. 

The phase diagram of the PAM is expected to have the form shown in Fig. 1. For D < D* 
the behavior along the critical line Ac(-D) should be that of DP, given that such behavior 
is generic for absorbing-state phase transitions without special symmetries or conserved 
quantities o], [l^. If < 1, then we expect mean-field like critical behavior as A — Ac = 
at fixed D > D*. Within the absorbing phase, for < A < Ac(-D), an isolated particle can 
produce an offspring, leading to annihilation of both the original and the new particle. On 
the line A = 0, this channel to annihilation is not available, and isolated particles cannot 
disappear. Thus the dynamics at long times, for D > 0, will be that of the diffusive 
annihilation process A + A ^ 0, for which the partic 
~ {In t)/t in two dimensions, and ~l/tin(i>3tll 



e density p(t) decays ~ 1/vt in 



12| . Finally, at the point X = D = 0, 



starting from all sites occupied, pairs are annihilated successively until only isolated particles 
remain. This is equivalent to the random sequential adsorption (RSA) of dimers. (In the 
present case, of course, dimers are removed not adsorbed, so the final particle density is 
1 — 26*00, where 6oo is the final coverage in RSA.) On the line, the final density of isolated 
particles is = 0.135335... 13|, while in two dimensions one has poo — 0.093108(8) fl^ . 



One may anticipate interesting crossover behaviors in the vicinity of one or another hmit. In 
the present work, however, we focus on determining the function Ac(-D) using Monte Carlo 
simulation. 




FIG. 1: Schematic phase diagram of the PAM in the X-D plane. The results of [71] imply that D* = 1 
in one and two dimensions. 



III. CLUSTER APPROXIMATIONS 



In this section we study the PAM through mean-field - site and pair approximations 15 |. 
In general, mean-field results provide a good qualitative description of the phase diagram 
and give an order-of-magnitude estimate of the critical point, n-site approximations are 
a natural way to improve the mean-field approach. The method consists of treating the 
transitions inside clusters of n sites exactly, while transitions involving sites outside the 
cluster are treated in an approximate manner. 



A. One-site approximation 



Let p = Prob((Ta; = 1) denote the density of particles. The density is governed by, 

e 



dp 

II- D 



- d5^P((t, = l,a,+e = 0) 

e 

+ Z}5^PK = 0,a,+, = 1), (2) 

e 

where the sums are over the 2d nearest-neighbors of site x, and P{(Jx, (Jx+e) is a two-site 
joint probabihty. Equation (|2]) couples the one-site probability p to the two-site proba- 
bilities, which in turn depend on the three-site probabilities, and so forth, leading to an 
infinite hierarchy of equations for the n-site probabilities. The site-approximation consists 
in truncating this hierarchy at n = 1, so that the two-site probabilities are replaced by a 
product of two one-site probabilities. Assuming spatial homogeneity and isotropy we obtain 
the following equation for p, 

|.i_^[A,_p,Ay] (3) 

The stationary solutions are p = (unstable for A > 0) and p = A/(2 + A). Thus, in 
this approximation the critical parameter Ac is zero in any dimension. Notice that in this 
approximation the diffusion rate has no influence on the stationary solution. 



B. Pair approximation 



To derive the pair approximation equations, we consider the changes in the configuration 
of a NN pair of sites (the central pair), given the states of the surrounding sites. Using the 
symbols o and • to represent, respectively, vacant and occupied sites, the states of a pair 
are oo, ••, wo, and o«; the latter two have the same probability and may be treated as a 
single class using appropriate symmetry factors. It is convenient to use (••) to denote the 
fraction of •• pairs, and so on. Then we have for the site fractions (•) = p and (o) = 1 — p: 



(•) = (••) + (»o), (4) 
(o) = (oo) + (.o). (5) 

The pair fractions satisfy (oo) + 2(o») + (««) = 1. The pair approximation consists in writing 
the joint probabihty of a set of three neighboring sites in the form (abc) = {ab){bc)/ (6). 

There are five possible transitions between the pair states. Consider for example the 
transition oo — > o«. This can occur via creation or via hopping, if and only if the rightmost 
site of the central pair has an occupied NN. Since its NN within the central pair is vacant, at 
least one of its 2d — 1 NNs outside the central pair must be occupied. The rate of transitions 
via creation is 

Ri,c - (1 - D)X^^ ^ (6) 

where we introduced A = A/(l + A). Adding the contribution due to diffusion, we obtain 
the total rate for this transition. 

Note that the contribution to the loss term for (oo) associated with this process is 2i?i, due 
to the mirror transition oo — > #0. 

The rates for the other transitions are: 

o» 00: 



o« — >■ ••: 



— >■ 00: 



^4 = ^(l-i^)(l-A)(..) (10) 



01 



2d - 1 (..) 



2d (•) 

The equations of motion for the pair probabilities are then 



2(1-L>)(1- A)(««) + D(o«) (11) 



-(00) = 2i?2 + i?4-2i?i 



(12) 



l(o«) =i?i+i?5-i?2-i?3 (13) 

and 

^(••) = 2i?3-i?4-2i?5 (14) 



The active stationary solution of the above equations is 

X[{4d - 3 + D)X - 2{1 - 2dD)] 

^ ~ (4rf - 3 + D)A2 + 2[2d{D + 2) - 3] A + A{2d - 1)D' ^ ^ 

and 

(16) 



For A < 2{l~2dD)/{4:d-3 + D), only the trivial solution (p = 0) exists, li D > D* = l/2d, 
however, the active solution exists for any A > 0. The phase transition occurs at 



'' 2{i~2dD) n <- n* = J- 

4d-3+D ' 1^ 1^ 2d 



Ac 



(17) 

0, D> D* 



Thus the pair approximation predicts a nonzero critical creation rate only for diffusion 
probabilities D < D* = l/{2d)] for larger values of D, there is a nonzero particle density 
for any A > 0, as in the one-site approximation. For D = 0, we have Ac = 2, 2/5, and 2/9 
in one, two and three dimensions, respectively; the corresponding values from simulation 
are Ac = 5.368(1), 1.0156(1), and 0.475(1). [We note that the pair approximation results 
derived above differ slightly from those given in ^\ since in the latter case the annihilation 
rate for a NN particle pair is taken as (1 — i^)/(l + A), i.e., d times the rate given in Eq. 

a 

Katori and Konno [T] proved that the prediction of D* < 1, furnished by the pair approxi- 
mation, is wrong for d < 2. That is, in one and two dimensions, Ac > for any D < 1. In the 
following section we investigate how Ac tends to zero as D — )■ 1 in one and two dimensions, 
and determine D* in the three-dimensional case. 



IV. SIMULATIONS 



We use Monte Carlo simulations to obtain accurate values of the critical creation rate 
Ac(-D) and the critical exponents of the PAM in one, two, and three dimensions. 

A. One dimension 

1. Spreading behavior 

A well established method for determining the critical point and certain critical exponents 
is through the study of propagation of activity, starting from a localized seed, as proposed 
long ago by Grassberger and de la Torre [l6|. One studies the activity in a large set of 
trials, all starting from a configuration very close to the absorbing state. Here the initial 
configuration is that of a single pair of particles at the two central sites, in an otherwise 
empty lattice. Each trial ends when it reaches the absorbing state, or at a maximum time, 
t^axi chosen such that the activity never reaches the edges of the system (in any trial) for 

For A > Ac there is a nonzero probability that the process survives as t — )► oo; for A < Ac 
the process dies with probability 1. Of primary interest are P{t), the probability of surviving 
until time t or greater, n{t), the mean number of particles (averaged over all trials), and 
R^{t), the mean-square distance of particles from the origin. At the critical point these 
quantities follow asymptotic power laws. 

Pit) oc (18) 
nit) oc t"^ (19) 
R^{t) oc e^^. (20) 

The exponents 5, r], and Zsp satisfy the hyperscaling relation 45 + 2?7 = dzsp., in c? < 4 



dimensions 16|. (We note that Zgp is related to the usual dynamic exponent z via Zgp = 2/z.) 

We study activity spreading in one dimension using samples of from 10^ or 2 x 10® trials 
for each A value of interest. The maximum time tmax = 15 000 for D < 0.7, 30 000 for 
D = 0.8 and 0.9, and 60 000 for D = 0.95. (As D increases, the asymptotic power-law 
behavior occurs at ever later times.) To ensure that activity never reaches the borders, we 
use a lattice size of L = 50 000 for tmax = 15 000, and L = 80 000 for the longer studies. 



A study performed at a given value of A is used to generate results for nearby values using 
sample reweighting [17]. 

To locate the critical point, we use the criterion of power-law behavior of n{t); Fig. [2] 
illustrates the analysis for D = 0.3. The main graph is a log-log plot of n{t) showing an 
apparent power law for A = 3.4687. The curves for nearby values (specifically, A = 3.4681, 
3.4684, 3.4690, and 3.4693, obtained via reweighting), cannot be distinguished on the scale 
of this graph, but if we plot n* = n/t^, the curves for different A values fan out (upper 
inset), with upward curvature indicating a supercritical value of A and vice- versa. 

The exponent 1] is estimated via analysis of the local slope, //(t), defined as the inclination 
of a least-square linear fit to the data (on logarithmic scales), on the interval [t/a, at]. (The 
choice of the factor a represents a compromise between high resolution, for smaller a, and 
insensitivity to fluctuations, for larger values; here we use a = 2.59.) Plotting rjit) versus 
1/t (lower inset of Fig. [2]) allows one to estimate Ac (the curves for A > Ac veer upward, 
and vice- versa), and to estimate the critical exponent rj by extrapolating r]{t) to 1/t —t- 0. 
The main source of uncertainty in the exponent estimates is the uncertainty in Ac itself. 
An analogous procedure is used to estimate exponents 6 and Zsp. In Table |T] we list the 
critical parameters and spreading exponents found via spreading simulations combined with 
local-slopes analysis. 

For D = 0, the critical parameter for the PAM, Ac(0)=5. 368(1), is considerably larger 
than that of the contact process (Ac=3. 29785(2)), as expected since here each annihilation 
event removes two particles. (The fact that Ac is less than twice the corresponding value in 
the CP may be attributed to the tendency for particles to cluster: removing two particles 
may eliminate additional pairs, thus reducing the effective rate of annihilation.) 

For all diffusion probabilities studied, our estimates for the critical exponents are in good 
accord with the DP values 6 = 0.15947(3), r] = 0.31368(4), and z = 1.26523(3) [2]. A plot 
of the phase boundary in the X-D plane (see Fig. [3]) suggests that Ac — )■ as -D — t- 1, so 
that D* = 1 in agreement with the Katori-Konno theorem. Extrapolation of D versus Ac, 
using a fourth-order polynomial fit to the data for D > 0.6, yields D = 1.0005 for Ac = 0, 
confirming to high precision that Ac > for D < 1. 




FIG. 2: Main graph: n{t) on log scales for the one-dimensional PAM with D = 0.3 and A = 3.4687. 
Upper inset: n* = n/f on log scales, for (lower to upper) A = 3.4681, 3.4684, 3.4687, 3.4690, and 
3.4693. Lower inset: local slopes rj{t) for the same set of A values. 

B. Two dimensions 

1. Steady- state behavior 

We encounter rather large uncertainties in studies of spreading behavior of the two- 
dimensional PAM, and so turn to the steady-state approach to investigate this system (on 
the square lattice). In these studies we initialize the system with all sites occupied, and 
allow it to evolve until it attains a quasistationary (QS) regime, in which bulk properties 
such as the particle density p, averaged over surviving realizations, are time- independent. 



According to the finite-size scaling hypothesis [l8|, [l9|, the QS properties depend on system 
size L through the ratio L/^^ or equivalently through the scaling variable AL^/''^, where 



D 


Ac 


5 


V 




0.0 


5.3720(5) 


0.159(1) 


0.315(1) 


1.266(2) 


0.1 


4.6709(2) 


0.161(1) 


0.314(2) 


1.264(1) 


0.2 


4.0417(2) 


0.160(1) 


0.314(1) 


1.266(1) 


0.3 


3.4687(2) 


0.162(1) 


0.312(1) 


1.268(2) 


0.4 


2.9411(2) 


0.160(2) 


0.314(2) 


1.264(2) 


0.5 


2.4473(1) 


0.159(1) 


0.315(2) 


1.267(3) 


0.6 


1.9778(2) 


0.159(1) 


0.315(2) 


1.266(1) 


0.7 


1.5231(2) 


0.158(2) 


0.315(1) 


1.265(2) 


0.8 


1.0684(2) 


0.159(2) 


0.315(2) 


1.265(2) 


0.9 


0.5891(1) 


0.161(1) 


0.315(1) 


1.267(3) 


0.95 


0.3214(1) 


0.159(2) 


0.318(3) 


1.266(4) 



TABLE I: Results of spreading simulations for the RAM in one dimension. 

A = A — Ac- Expressing the order parameter as a function of A and L, we have 

p(A, L) oc L-^'^^fi^L^'^^). (21) 
with f\x) oc as X — )■ oo. At the critical point, A = 0, 

p(0,L) oc L-^/*^^. (22) 

Thus an asymptotic power-law dependence of p on L is a useful criterion for criticality. 

We study the QS density as a function of system size to locate the critical point, using 
sizes L = 25, 50, 100,..., 800. The relaxation time varies from r = 800 for the smallest size, to 
r = 200 000 for the largest; the number of realizations varies from 500 to 10 000. Using the 
power-law criterion, we obtain the estimates for Ac listed in Table [TTl It is worth mentioning 
that the values for Ac, for D = and D = 0.1 are in good agreement with those obtained in 
the preliminary spreading behavior studies. 

In Fig. m we verify the scaling collapse of the order parameter, plotting x = L^^'^^p versus 
y = AL^/*^^, for system sizes L=16, 32, 64, 128, and 256. A good collapse is obtained using 
the DP values = 0.733 and P/v^ = 0.795 12||. The data are consistent with the scaling 
law p oc A'^, using the DP value /3 = 0.583(4) |2|. 
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FIG. 3: Points along the critical line in the X-D plane in one dimension, as determined via simulation. 
Error bars are smaller than symbols. The solid line is a quartic fit to the six points with largest D. 



D 


Ac 


0.0 


1.0156(1) 


0.1 


0.7877(1) 


0.2 


0.5890(5) 


0.3 


0.4166(1) 


0.4 


0.2685(5) 


0.5 


0.1462(2) 


0.6 


0.056(1) 



TABLE II: Critical parameters obtained through steady-state simulations in two dimensions. 

2. Quasistationary simulations 

As D approaches 0.7 the critical value Ac becomes very small. We require an efficient 
simulation method to obtain precise estimates for the critical value for larger diffusion rates, 
in particular, to determine how Ac tends to zero as D increases. For this purpose we use 
quasistationary (QS) simulations, which sample directly the QS probability distribution, 
that is, the long-time distribution conditioned on survival. The details of the method are 



In X 



FIG. 4: (Color online) Scaling plot of the stationary density in the two-dimensional PAM with D = 0. 
System sizes L=16 (+); 32 (x); 64 (diamonds); 128 (•) and 256 (squares). The slope of the straight 
line is 0.583. 



explained in Ref. 



20|]. To obtain these results we use lattice sizes L = 100, 200, 400 and 



800 for D = 0.7, and include studies of larger systems for higher diffusion rates (up to 
L = 6400, for D > 0.78). The critical point is determined via the criteria of power-law 
scaling of the density and mean lifetime with system size, and convergence of the moment 
ratio m = (p^)/p^ to a finite limiting value as L — oo, as discussed in 2l(|. (The lifetime r 
is expected to follow r ~ L^.) Using this method we obtain the values listed in Table UTTl 
We note that our results for z, and the limiting moment ratio rric are consistent with 



13, 



2l|. 



the known DP values of 0.795(10), 1.7674(6), and 1.3257(5), respectively 
For D = 0.8, Ac is of order 10~^, and a precise determination becomes very difficult due to 
the small number of particles present in the system. Reliable determination of Ac for larger 
diffusion rates would therefore require studies of even larger systems, which was deemed 
impractical. 

We find that Ac(-D) can be fit quite well using an expression of the form 



D 




0.60 


0.05632(3) 


0.70 


0.00940(5) 


0.73 


0.003957(3) 


0.78 


0.0004815(7) 


0.80 


0.00015(2) 



TABLE III: Critical parameters obtained through quasistationary simulations in two dimensions. 



Ac = Aexp 



C 



(23) 



Applied to the data for D > 0.4, a least-squares procedure yields 7 = 1.41(2), C = 0.984(2), 
and A = 2.02(2). The good quality of the fit is evident in the inset of Fig. [5l Thus, while a 
plot of the data on linear scale might suggest that Ac — J- at some diffusion rate between 0.8 
and 1 (see Fig. [5l main graph), our results are in fact consistent with Ac nonzero, though 
very small, for diffusion rates between 0.7 and unity. 



C. Three dimensions 



We employed quasistationary simulations to determine Ac(-D) for the PAM on the simple 
cubic lattice. For relatively small diffusion rates good results are obtained using lattice sizes 
L = 8, 16, 24, 36, and 54. For diffusion rates greater than about 0.25, however, there are 
substantial finite-size effects, and to observe clear signs of DP-like scaling we need to study 
larger systems {L = 80 and 120 in addition to the sizes mentioned above). The results (see 
Table HVl and Fig. |6]), show that in this case Ac does fall to zero at a diffusion rate considerably 
less than unity; extrapolation of the data to A = yields D* = 0.333(3). The critical 
exponents determined via finite-size scaling analysis, P/iy^ = 1.40(1) and z = 1.94(2), are 
once again in good agreement with hterature values of 1.39(3) and 1.919(4), respectively. Our 
study yields the moment ratio value m = 1.47(1) for the three-dimensional models in the DP 
universality class; to our knowledge this quantity has not been determined previously. For 
D > D*, the particle density is expected to tend to zero linearly with A, as the reproduction 
rate approaches zero. We have verified this behavior (down to A = 10~^) for D = 0.8. 




X 



FIG. 5: Critical line of the two-dimensional PAM. Inset: the same data plotted as In Ac versus 

1/(1 



D 


Ac 


0.0 


0.47390(5) 


0.1 


0.2943(1) 


0.2 


0.1420(1) 


0.25 


0.07790(5) 


0.28 


0.04487(3) 


0.31 


0.01762(2) 


0.32 


0.0103(1) 



TABLE IV: Critical parameters obtained through quasistationary simulations in three dimensions. 
V. DISCUSSION 



We study the phase boundary of the pair annihilation model in the reproduction rate - 
diffusion probability (A - D) plane. Our simulation results are consistent with the theorem 



0.0 0.1 0.2 0.3 0.4 0.5 



FIG. 6: Critical line of the PAM in three dimensions; error bars are smaller than symbols. The solid 
line is a cubic fit to the data, yielding D* = 0.333(3). 

proven some time ago by Katori and Konno [7], namely that in one and two dimensions, 
Ac > for any D < 1. The pair approximation is in conflict with this result, as it predicts 
that in any number of dimensions, there is a diffusion probability D* < 1, above which 
Ac = 0. In one dimension the behavior (in simulations) is straightforward, as Ac oc 1 — D 
for D ~ 1. In two dimensions however it is quite subtle, as Ac becomes exponentially small 
as D — )■ 1, and a cursory analysis could well give the impression that Ac is actually zero at 
some value of D between 0.8 and unity. Finally in three dimensions the pair approximation 
prediction is verified qualitatively; we find D* = 0.333(3) in this case, while the PA yields 
D* = 1/6. Intuitively, the unusual behavior of Ac(-D) in two dimensions can be understood 
as a consequence of d = 2 marking a critical dimension for the recurrence of a random 
walk. Our simulation results for critical exponents and the moment ratio m are consistent 
with the directed percolation values, as expected. Given the qualitative failure of the pair 
approximation in one and two dimensions, it is natural to ask whether approximations using 
larger clusters would predict the phase diagram correctly. This strikes us as unlikely, since 



cluster approximations have been found to be insensitive to subtle effects involving diffusion 
and/or multiparticle rules in other cases 15|, |22|, |23 |. 
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